Covid19 Data in Japan
日本の新型コロナウィルスに関する情報が抱える大きな問題の一つに、ソースごとの値が微妙にことなっていることが挙げられる。原因は明白で、総元締めであるはずの厚生労働省が各自治からの患者データを集計せずに各自治体のプレスリリースやウェブサイトの数値を積上げたものを発表しているからである。
本来であれば医療機関->市町村->都道府県->厚生労働省が一気通貫して同じデータを扱えるようにすべきであり、このあたりの改革は、河野行政改革担当大臣に期待するところである。
さて、理想論はさておき、実際のデータはどのようになっていて、そこから何が読み取れるかを見てみる。
厚生労働省データ
厚生労働省のオープンデータのページで公開しているデータは、あくまでも、「国内事例」のみであり、チャーター便、空港検疫、外国船籍の客船を除くデータになっている。なお、データによっては途中からデータソースが変更されているものがあるため可視化すると不自然な変動が見られる場合がある。
Import
項目ごとに個別のCSVファイルになっている。
陽性者数
PCR検査で陽性が確定した「陽性者数」のデータ。厚生労働省が公開している陽者数のデータは単日のデータのみ。
"https://www.mhlw.go.jp/content/pcr_positive_daily.csv" %>%
readr::read_csv()
死亡者数
死亡者数のデータは、なぜか累計データだけが公開されている。
"https://www.mhlw.go.jp/content/death_total.csv" %>%
readr::read_csv()
入院治療等を要する者の数
いわゆる無症状患者などの軽症患者を除く陽性者の数。各報告日時点の集計値となっている。
"https://www.mhlw.go.jp/content/cases_total.csv" %>%
readr::read_csv()
PCR検査実施人数
件数ベースと人数ベースのデータが混在しているため数値が過大な傾向になっている。陽性者数と同じく単日。ただし、異なる数を合算しているので目安程度にしかならないと思われる。
"https://www.mhlw.go.jp/content/pcr_tested_daily.csv" %>%
readr::read_csv()
Wrangle (tidy and transform)
各データを扱いやすいように変換する。
陽性者数
単日のデータなので累計データを求める。なお、記録漏れがあっても対処できるように日付の補完を行う。
df_mhlw <- "https://www.mhlw.go.jp/content/pcr_positive_daily.csv" %>%
readr::read_csv() %>%
dplyr::mutate(date = lubridate::as_date(`日付`)) %>%
tidyr::complete(date = seq.Date(from = min(date), to = max(date), by = "day")) %>%
dplyr::mutate(n = as.integer(`PCR 検査陽性者数(単日)`),
cum = cumsum(n)) %>%
dplyr::select(date, n, cum)
df_mhlw
死亡者数
陽性者数とは逆に短日の数値を求める。
dead_mhlw <- "https://www.mhlw.go.jp/content/death_total.csv" %>%
readr::read_csv() %>%
dplyr::mutate(date = lubridate::as_date(`日付`)) %>%
tidyr::complete(date = seq.Date(from = min(date), to = max(date), by = "day")) %>%
dplyr::mutate(cum_dead = as.integer(`死亡者数`),
dead = cum_dead - dplyr::lag(cum_dead, default = 0)) %>%
dplyr::select(date, dead, cum_dead)
dead_mhlw
入院治療等を要する者の数
各報告日時点の集計値。
patient_mhlw <- "https://www.mhlw.go.jp/content/cases_total.csv" %>%
readr::read_csv() %>%
dplyr::mutate(date = lubridate::as_date(`日付`)) %>%
tidyr::complete(date = seq.Date(from = min(date), to = max(date), by = "day")) %>%
dplyr::mutate(patient = as.integer(`入院治療を要する者`)) %>%
dplyr::select(date, patient)
patient_mhlw
PCR検査実施人数
pcr_mhlw <- "https://www.mhlw.go.jp/content/pcr_tested_daily.csv" %>%
readr::read_csv() %>%
dplyr::mutate(date = lubridate::as_date(`日付`),
pcr = as.integer(`PCR 検査実施件数(単日)`)) %>%
dplyr::select(date, pcr) %>%
tidyr::complete(date = seq.Date(from = min(date), to = max(date), by = "day"),
fill = list(pcr = 0L)) %>%
dplyr::mutate(cum_pcr = cumsum(pcr))
pcr_mhlw
結合
上記のデータフレームをひとつにまとめる。
x_mhlw <- df_mhlw %>%
dplyr::left_join(dead_mhlw, by = c("date")) %>%
dplyr::left_join(patient_mhlw, by = c("date")) %>%
dplyr::left_join(pcr_mhlw, by = c("date")) %>%
tidyr::replace_na(replace = list(dead = 0L, cum_dead = 0L,
patient = 0L,
pcr = 0L, cum_pcr = 0L)) %>%
dplyr::mutate(diff = n - dplyr::lag(n, default = 0L),
diff_dead = dead - dplyr::lag(dead, default = 0L),
diff_patient = patient - dplyr::lag(patient, default = 0L),
diff_pcr = pcr - dplyr::lag(pcr, default = 0L)) %>%
dplyr::select(date, n, diff, cum,
patient, diff_patient,
dead, diff_dead, cum_dead, pcr, diff_pcr, cum_pcr)
x_mhlw
Visualize
陽性者数の推移
最近ではPCR検査の件数が大きく増えており、その感度が一説には約35%と言われていることからみても、陽性者数だけで議論することは意味がないと言われているが、まずは可視化する。
scale_axis <- 50
x_mhlw %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity",
width = 1.0, alpha = 0.75) +
ggplot2::geom_line(ggplot2::aes(y = cum / scale_axis), colour = "dark green") +
ggplot2::scale_y_continuous(
name = "陽性者数（単日）",
sec.axis = ggplot2::sec_axis(~ . * scale_axis,
name = "陽性者数累計（折線）")
)
2020年10月29日の報道では陽性者数の累計が10万人を超えた（チャーター便、空港検疫など含む）と報道されていたが、2020-11-02 21:55:45 時点で厚生労働省が公開しているデータに基づくと 2020-11-01 の \(101368\) 人である。このラグは集計方法に起因すると推定できる。
死亡者数の推移
陽性者数の推移に変わり着目されているのが死亡者数の推移である。まずは、単純な死亡者数の推移。
scale_axis <- 20
x_mhlw %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::geom_bar(ggplot2::aes(y = dead), stat = "identity",
width = 1.0, alpha = 0.75) +
ggplot2::geom_line(ggplot2::aes(y = cum_dead / scale_axis),
colour = "dark green") +
ggplot2::scale_y_continuous(
name = "死者数（単日）",
sec.axis = ggplot2::sec_axis(~ . * scale_axis,
name = "死者数累計（折線）")
)
plotly::ggplotly()
入院治療等を要する者の推移
x_mhlw %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::geom_bar(ggplot2::aes(y = patient), stat = "identity",
width = 1.0, alpha = 0.75) +
ggplot2::labs(y = "要入院治療者数（集計）")
値が急激に変わっている部分が認められるので、前日比のデータをプロットしてみる。
gg_pcr <- x_mhlw %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::geom_line(ggplot2::aes(y = diff_patient), colour = "dark green") +
ggplot2::labs(y = "要入院治療者数（前日比）")
gg_pcr %>%
plotly::ggplotly()5月8日の時点で急激に値が変わっているのが分かる。これは、オープンデータのページに記載されているように「令和2年（2020年）5月8日公表分から、データソースを従来の厚生労働省が把握した個票を積み上げたものから、各自治体がウェブサイトで公表している数等を積み上げたものに変更した。」ためである。
PCR検査実施人数の推移
scale_axis <- 20
x_mhlw %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::geom_bar(ggplot2::aes(y = pcr), stat = "identity",
width = 1.0, alpha = 0.75) +
ggplot2::geom_line(ggplot2::aes(y = cum_pcr / scale_axis),
colour = "dark green") +
ggplot2::scale_y_continuous(
name = "PCR検査実施人数（単日）",
sec.axis = ggplot2::sec_axis(~ . * scale_axis,
name = "累積人数（折線）")
)
陽性者数 vs 検査人数
検査人数が増えたから陽性者数も増えたと言われているが、本当なのか確認のためのプロットしてみる。
scale_axis <- 100
x_mhlw %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::geom_line(ggplot2::aes(y = n), colour = "dark green") +
ggplot2::geom_line(ggplot2::aes(y = pcr / scale_axis), color = "dark blue") +
ggplot2::scale_y_continuous(
name = "陽性者数（単日）",
sec.axis = ggplot2::sec_axis(~ . * scale_axis,
name = "検査人数（単日）")
)
「令和2年9月29日より、福岡県における民間検査機関等での実施件数（過去行われたものも含む）を計上している。」ため異常値になっている9月29日のPCR検査実施人数のデータを削除しておく。そもそも、過去のデータを単日で積上げてしまうやり方は非常に疑問である。
傾向を把握しやすくするために両者とも7日間の移動平均も合わせてプロットする。
scale_axis <- 10
x_mhlw %>%
dplyr::filter(date != lubridate::as_date("2020-09-29")) %>%
tidyr::complete(date = seq.Date(from = min(date), to = max(date), by = "day")) %>%
dplyr::mutate(ma = zoo::rollmeanr(n, k = 7L, na.pad = TRUE),
pcr_ma = zoo::rollmeanr(pcr, k = 7L, na.pad = TRUE)) %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::geom_line(ggplot2::aes(y = n), colour = "light green") +
ggplot2::geom_line(ggplot2::aes(y = pcr / scale_axis), color = "light blue") +
ggplot2::geom_line(ggplot2::aes(y = ma), colour = "dark green") +
ggplot2::geom_line(ggplot2::aes(y = pcr_ma / scale_axis), color = "dark blue") +
ggplot2::scale_y_continuous(
name = "陽性者数（単日）濃緑：移動平均",
sec.axis = ggplot2::sec_axis(~ . * scale_axis,
name = "検査人数（単日）濃青：移動平均")
)## Warning: Removed 6 row(s) containing missing values (geom_path).
## Warning: Removed 6 row(s) containing missing values (geom_path).
このグラフから分かるように7月〜8月の検査人数の増加に伴い陽性者数も増加、その後の減少・横ばいに伴い陽性者数も減少・横ばいしていることが読み取れる。
陽性者数 vs 要入院治療者数
scale_axis <- 0.2
x_mhlw %>%
dplyr::filter(date != lubridate::as_date("2020-04-22")) %>%
tidyr::complete(date = seq.Date(from = min(date), to = max(date), by = "day")) %>%
dplyr::mutate(ma = zoo::rollmeanr(cum, k = 7L, na.pad = TRUE),
patient_ma = zoo::rollmeanr(patient, k = 7L, na.pad = TRUE)) %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::geom_line(ggplot2::aes(y = cum), colour = "light green") +
ggplot2::geom_line(ggplot2::aes(y = patient / scale_axis), color = "light blue") +
ggplot2::geom_line(ggplot2::aes(y = ma), colour = "dark green") +
ggplot2::geom_line(ggplot2::aes(y = patient_ma / scale_axis), color = "dark blue") +
ggplot2::scale_y_continuous(
name = "陽性者数（累計）濃緑：移動平均",
sec.axis = ggplot2::sec_axis(~ . * scale_axis,
name = "要入院治療者数（集計）濃青：移動平均")
)## Warning: Removed 6 row(s) containing missing values (geom_path).
## Warning: Removed 6 row(s) containing missing values (geom_path).
plotly::ggplotly()
陽性者数の急増に伴い要入院治療者数も杖ているが、9月以降の要入院治療者数は増えておらず10月以降は横ばい傾向にあることが分かる。つまり、陽性と判断されても入院治療を必要としない軽症もしくは無症状の割合が増えていると考えられる。
陽性者数 vs 死亡者数
陽性者数の増加は検査人数の増加にほぼ比例していると見ることができるが、死者数のはどうであろうか？4月22日は異常値として除いておく。
scale_axis <- 0.1
x_mhlw %>%
dplyr::filter(date != lubridate::as_date("2020-04-22")) %>%
tidyr::complete(date = seq.Date(from = min(date), to = max(date), by = "day")) %>%
dplyr::mutate(ma = zoo::rollmeanr(n, k = 7L, na.pad = TRUE),
dead_ma = zoo::rollmeanr(dead, k = 7L, na.pad = TRUE)) %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::geom_line(ggplot2::aes(y = n), colour = "light green") +
ggplot2::geom_line(ggplot2::aes(y = dead / scale_axis), color = "light blue") +
ggplot2::geom_line(ggplot2::aes(y = ma), colour = "dark green") +
ggplot2::geom_line(ggplot2::aes(y = dead_ma / scale_axis), color = "dark blue") +
ggplot2::scale_y_continuous(
name = "陽性者数（単日）濃緑：移動平均",
sec.axis = ggplot2::sec_axis(~ . * scale_axis,
name = "死者数（単日）濃青：移動平均")
)## Warning: Removed 6 row(s) containing missing values (geom_path).
## Warning: Removed 6 row(s) containing missing values (geom_path).
plotly::ggplotly()
陽性者が増加してから3〜4週間ほど遅れて死亡者数が増える傾向が、4月と8月の陽性者数のピークからよも取れる。しかし、両ピークを比較すると8月のピーク時の方が陽性者数に対する死者数の割合が低くなっていると言える。
要入院治療者数 vs 死亡者数
scale_axis <- 0.001
x_mhlw %>%
dplyr::filter(date != lubridate::as_date("2020-04-22")) %>%
tidyr::complete(date = seq.Date(from = min(date), to = max(date), by = "day")) %>%
dplyr::mutate(patient_ma = zoo::rollmeanr(patient, k = 7L, na.pad = TRUE),
dead_ma = zoo::rollmeanr(dead, k = 7L, na.pad = TRUE)) %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::geom_line(ggplot2::aes(y = patient), colour = "light green") +
ggplot2::geom_line(ggplot2::aes(y = dead / scale_axis), color = "light blue") +
ggplot2::geom_line(ggplot2::aes(y = patient_ma), colour = "dark green") +
ggplot2::geom_line(ggplot2::aes(y = dead_ma / scale_axis), color = "dark blue") +
ggplot2::scale_y_continuous(
name = "要入院治療者（集計）濃緑：移動平均",
sec.axis = ggplot2::sec_axis(~ . * scale_axis,
name = "死者数（単日）濃青：移動平均")
)## Warning: Removed 6 row(s) containing missing values (geom_path).
## Warning: Removed 6 row(s) containing missing values (geom_path).
plotly::ggplotly()
当然と言えるが、要入院治療者が増えると、その数週間後に死亡者数が増えていることが分かる。